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Abstract 

The Last Interglacial (LIG) stage (ca. 130-115 ka), with polar temperatures likely 3-5° C warmer than today, 
serves as a partial analogue for low-end future warming scenarios. Multiple indicators suggest that LIG global sea 
level (GSL) was higher than at present; based upon a small set of local sea level indicators, the Intergovernmental 
Panel on Climate Change (IPCC)'s Fourth Assessment Report inferred an elevation of approximately 4-6 m. While 
this estimate may be correct, it is based upon overly simplistic assumptions about the relationship between local 
sea level and global sea level. Sea level is often viewed as a simple function of changing global ice volume. This 
perspective neglects local variability, which arises from several factors, including the distortion of the geoid and the 
elastic and isostatic deformation of the solid Earth by shifting ice masses. Accurate reconstruction of past global 
and local sea levels, as well as ice sheet volumes, therefore requires integrating globally distributed data sets of local 
sea level indicators. To assess the robustness of the IPCC's global estimate and search for patterns in local sea 
level that are diagnostic of meltwater sources, we have compiled a comprehensive database that includes a variety 
of local sea level indicators from 47 localities, as well as a global sea level record derived from oxygen isotopes. 
We generate a global synthesis from these data using a novel statistical approach that couples Gaussian process 
regression to Markov Chain Monte Carlo simulation of geochronological errors. Our analysis strongly supports the 
hypothesis that global sea level during the Last Interglacial was higher than today, probably peaking between 6-9 
m above the present level. This level is close to that expected from the complete melting of the Greenland Ice 
Sheet, or from major melting of both the Greenland and West Antarctic Ice Sheets. In the period when sea level 
was within 10 m of the modern value, the fastest rate of sea level rise sustained for a 1 ky period was likely about 
80-110 cm per century. Combined with the evidence for mildly higher temperatures during the LIG, our results 
highlight the vulnerability of ice sheets to even relatively low levels of sustained global warming. 
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1. Introduction 

As a result of industrial activity, greenhouse gas 
concentrations now exceed levels reached on Earth at 
any time within the last eight hundred thousand years 
(Solomon et al., 2007). The long-lived greenhouse gases 
in the atmosphere today produce a radiative forcing 
of about 2.6 W/m^. Given a climate sensitivity of 2- 
4.5°C per doubling of carbon dioxide levels (Meehl et al., 
2007), the equilibrium global warming expected from 
this forcing - without considering the effects of any fur- 
ther increases in greenhouse gas concentrations - is be- 
tween 1.4-3.2°C. Among the many effects expected to 
accompany this warming, global sea level rise, driven 
primarily by thermal expansion of seawater and melt- 
ing ice sheets, features prominently (Meehl et al., 2007). 
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Uncertainties in ice sheet behavior make it difficult to 
predict sea level rise precisely, but by the end of the 
twenty-first century, mean global sea level rise could ex- 
ceed one meter (e.g., Rahmstorf, 2007; Grinsted et al., 
2009). Since changes of this magnitude have no prece- 
dent in human experience, to understand them and to 
compile observations against which to test models of 
future climate change, it is necessary to turn to the ge- 
ological record. 

Data from past interglacial stages and earlier warm 
periods may help constrain the future behavior of sea 
level. Synthesizing geological sea level indicators into a 
global reconstruction requires taking into account inter- 
pretative uncertainties, geochronological uncertainties, 
and regional variability. Differences between local and 
global sea level arise from several factors, including the 
distortion of the geoid and the elastic and isostatic de- 
formation of the solid Earth by shifting ice masses (Far- 
reh and Clark, 1976; Mitrovica and Milne, 2003). Ac- 
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curate reconstruction of past global and local sea levels 
therefore requires integrating global data sets of local 
sea level indicators. 

In this paper, we develop a novel statistical approach 
to assessing the probability distribution of local and 
global sea level through time from a database of geo- 
graphically distributed sea level indicators. Our algo- 
rithm is based upon Gaussian process regression with 
a Markov Chain Monte Carlo treatment of geochrono- 
logical errors. We apply this method to a new database 
of geographically dispersed sea level indicators spanning 
the Last Interglacial (LIG) stage, which climaxed about 
125 thousand years ago. The LIG (also known as the 
Eemian stage, its local northern European name, and as 
Marine Isotope Stage 5e) is of special interest for three 
reasons: (1) it is recent enough that it may be possible to 
construct relatively high-resolution sea level records for 
this interval, (2) due in large part to enhanced northern 
hemisphere insolation, mean global temperature may 
have been slightly warmer than at present, and (3) sev- 
eral lines of evidence suggest global sea level was higher 
than today, perhaps by 4-6 m (Jansen et al., 2007). 

Greenhouse gas concentrations during the LIG were 
comparable to pre-Industrial Holocene levels (Petit 
et al., 1999), but Earth's orbital eccentricity during the 
Last Interglacial was 0.038-0.041, more than twice as 
high as the modern value (Berger and Loutre, 1991). 
Energy balance modeling predicts that, as a conse- 
quence, summer temperatures between 132-124 ka on 
all land masses except Antarctica were at least 0.5° C 
warmer than today (Crowley and Kim, 1994), while a 
more complete climate model indicates summer temper- 
atures 2-4° C warmer than today in most of the Arc- 
tic (Otto-Bliesner et al., 2006). Ice core data from 
both Greenland and Antarctica suggest polar temper- 
ature anomalies in both hemispheres of about 3-5° C 
(Jansen et al., 2007), comparable to the 3-6° C of Arctic 
warming expected to accompany 1-2° C of global warm- 
ing (Katsov et al., 2004). In Europe, pollen data sug- 
gest middle Eemian summer temperatures about 2°C 
warmer than present (Kaspar et al., 2005). While the 
global mean temperature anomaly is uncertain, sea- 
surface temperatures in the equatorial Pacific (Lea, 
2004) and Atlantic (Weldeab et al., 2007) were about 
2°C warmer than pre-Industrial levels. 

For our analysis of LIG sea level, we compiled an ex- 
tensive database of sea level records from this time pe- 
riod. Sea level can be inferred from two basic sources: 
the oxygen isotopic composition of seawater as preserved 
in carbonates, and geomorphological and sediment olog- 
ical records of local sea level. The oxygen isotopic com- 
position of seawater is in large part a function of ice 
sheet volume and is thus closely related to global sea 
level. The latter is defined as the integrated local sea 
level over the entire ocean, which is a measure of ocean 
volume. Water enriched in the light isotope ^^O is con- 
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centrated preferentially in the ice sheets; therefore, as 
ice sheets melt and global sea level rises, the relative 
concentration of the heavy isotope, 5^^0, of seawater 
decreases. The record of marine oxygen isotopes is thus 
a proxy for global sea level, although its interpretation 
is complicated by several factors. The oxygen isotopic 
composition of marine water is also a function of local 
salinity, while the difference between the isotopic com- 
position of water and the carbonates precipitated from 
it is a function of temperature (Emiliani, 1966). The 
interpretation of the carbonate oxygen isotope record 
is further complicated by the histories of different wa- 
ter masses, by changes in their distribution over time, 
and, on long time scales, by isotopic exchange with the 
mantle (Jean-Baptiste et al., 1997; Schrag et al., 2002). 
Nonetheless, with assumptions about the isotopic com- 
position of ice sheets and a paleotemperature constraint, 
it is possible to estimate ice volume and thus sea level 
from S^^O in carbonate, as we do later. 

Many variables affect the local relative sea level 
recorded in sedimentary facies and in geomorphology. 
Aside from ocean volume, these factors include the ef- 
fect of the distribution of ice, water, and sediment on 
the geoid, lithospheric fiexure, isostatic adjustments, 
the orientation and magnitude of Earth's spin axis, and 
shoreline geometry (Farrell and Clark, 1976; Mitrovica 
and Milne, 2003), as well as tectonic uplift and ther- 
mal subsidence. As a result of gravitational and elastic 
effects, local sea level shortly after an ice sheet melt- 
ing event falls in the near-field and rises by more than 
the global average in the far-field. Over the follow- 
ing several millennia, isostasy slowly emplaces a mass 
of mantle to compensate for the missing ice, thereby 
changing the geoid so that local sea level in both the 
far-field and the near-field relax toward the global av- 
erage. (In the near-field, however, the vertical motion 
of the solid Earth complicates the observed relative sea 
level change.) For this reason, local sea levels at Pacific 
islands in the far-field of Laurentide Ice Sheet were 1-3 
m higher in the middle Holocene than today (Mitro- 
vica and Milne, 2002). Thus, even if LIG ice volume 
never shrunk below present levels and mean global sea 
level never exceeded its present value, local sea levels 
several meters higher than present could have occurred 
in the far-field of the Laurentide Ice Sheet early in the 
LIG, and comparably high local sea levels could have oc- 
curred in its intermediate-field late in the LIG (Lambeck 
and Nakada, 1992). Without accurate and precise dat- 
ing of the relevant sea level indicators, these effects could 
produce the false appearance of a magnified global sea 
level high-stand. In order to estimate ice sheet history 
from sea level records, it is thus necessary to account 
for these factors. Conversely, because many of these ef- 
fects lead relative sea-level changes to differ in the near- 
and far-fields of an ice sheet, a global database of local 
sea level indicators could address not just whether ice 
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volume was smaller during the Last Interglacial than 
today, but also what combination of melting ice sheets, 
if any, was responsible for higher global sea levels. 

Our goal is to use as comprehensive a database as 
possible to produce probability distributions of sea level 
as a function of geographic location and geological age. 
We have to cope with variable temporal uncertainty, 
as well as with variable errors in sea level measure- 
ments. In addition, some of the data are censored in 
that they provide only an upper or lower bound to sea 
level. Where possible, we also want to take advantage of 
quasi-continuous sequences, in which relative timing is 
known with greater precision than absolute dates. This 
is the case for the global oxygen isotope curve from ben- 
thic foraminifera (Lisiecki and Raymo, 2005), as well as 
for sequences of local observations of sedimentary facies 
from the Netherlands (Zagwijn, 1983) and of sea lev- 
els derived from hydrological modeling of foraminiferal 
oxygen isotopes from the Red Sea (Rohling et al., 2008). 

Our statistical approach is based upon Gaussian pro- 
cess regression (Rasmussen and Williams, 2006). This 
method, of which the commonly-used geospatial tech- 
nique of kriging interpolation is a well-known example, 
treats a field (such as sea level) as a collection of random 
variables drawn from a Gaussian distribution. By spec- 
ifying the covariance structure of the field, knowledge 
about the relevant physics affecting the process can be 
incorporated into the modeling without constraining it 
to fit a particular forward model. With a sufficiently 
precise and accurate data set, such an analysis will al- 
low us not only to place robust constraints on global 
sea level but also to identify the "fingerprints" produced 
by the gravitational, elastic, and isostatic effects of dif- 
ferent melting ice sheets (e.g., Mitrovica et al., 2001). 
It can thereby provide an independent test for different 
melt water sources in the Last Interglacial, and by ex- 
tension the possible susceptibility of each ice sheet to 
future melting. 

2. Database of LIG Sea Level Indicators 

We characterize each sea level indicator in our 
database by five parameters: its geographical position, 
its altitude with respect to mean tide level, its age, the 
range of depths at which it might have formed, and the 
local uplift or subsidence rate. With the exception of 
geographical position, each of these variables has uncer- 
tainties that we assume follow a Gaussian distribution. 
For some values, including all depositional depth ranges, 
uniform distributions between two limits a and b may be 
a better choice than Gaussian ones. In these cases, we 
substitute a Gaussian distribution with the same mean 
and standard deviation as the uniform distribution, i.e. 
(6-a)/Vl2. 

The full database is available on request from the au- 
thors and will be published with the final version of this 
paper. 
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2.1. Nature of the indicators and depositional ranges 

The sea level indicators take a variety of forms, in- 
cluding: constructional coral terraces that provide both 
geomorphological and ecological information; coral bio- 
facies in limestones that provide ecological but not ge- 
omorphological information; erosional features such as 
wave-cut terraces, sea caves, bioerosional notches, and 
raised beaches; and sediment ological and biofacial indi- 
cators of depositional depth. 

Most of the indicators refiect deposition or formation 
within a specific range of depths. The most common 
reef terraces and associated coral assemblages, for in- 
stance, are generally interpreted as indicating deposi- 
tion between mean low tide level and 5 m below mean 
low tide level (Lighty et al., 1982; Camoin et al., 2001). 
Intertidal sedimentary facies indicate deposition within 
the tidal range. While recognizing that LIG tidal am- 
plitudes could have been slightly different than today, 
we convert descriptive ranges such as these into a com- 
mon reference frame based on the tidal ranges reported 
in tide tables at a nearby modern locality. We also 
attempt to correct for variability in the measurement 
datum; while most sea level indicators have altitudes 
reported with respect to "modern sea level", some are 
more usefully described with reference to datums such 
as the mean low tide level or mean high tide level. We 
convert such datums into a mean tide level datum. 

Some data, such as subtidal sedimentary facies, are 
limiting points; they place an upper or lower limit on 
past sea level but do not indicate a specific depositional 
depth. In statistical terminology, limiting points are 
censored data. 

2.2. Age 

Age constraints on our data come from a variety of 
sources with a range of precisions. In some cases, age 
is constrained only by stratigraphic relationships with 
other units. In many cases, particularly involving coral 
reefs, radiometric (U/Th) dates are available. Other age 
constraints are derived from amino acid racemization, 
electron spin resonance dating, and related techniques 
such as thermoluminescence. 

In three cases (the global oxygen isotope curve, the 
Red Sea oxygen isotope curve, and the Dutch sea level 
curve) , relative ages are known with more precision than 
absolute ones. As described in the Appendix, we have 
scaled and shifted the age models of the Red Sea and 
Dutch local sea level curves to be consistent with the 
Lisiecki and Raymo (2005) age model for the global oxy- 
gen isotope curve. All of the dates outputted by our 
analysis should therefore be viewed within the context 
of this age model, which places the start of the Penul- 
timate Termination at 135 ka and the peak of the Last 
Interglacial at about 122-126 ka. 

When only a single conventional U/Th measurement 
from a unit is available, we expand the quoted ranges 
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Table 1: Sites, Number, and Types of Sea Level Indicators in the LIG Database 



Site 


# Obs. 


Type 


Reference 


Global 








Global oxygen isotope stack 


51 


isotopic 


Lisiecki and Raymo (2005) 


Northeastern Atlantic Ocean and Mediterranean 


ea 




Southern England 


2 


erosional 


Westaway et al. (2006) 


Bristol Channel, Britain 


1 


erosional 


Allen (2002) 


Belle Hogue Cave, Jersey 


1 


erosinal 


Keen et al. (1981) 


Port-Racme Beach, France 


1 


erosional 


Bates et al. (2003) 


The Netherlands 


8 


facies 


Zagwijn (1983) 


Hergla South, Tunisia 


2 


facies 


Hearty et al. (2007) 


Quaternary Basin, Mauretania 


2 


facies 


Giresse et al. (2000) 


Northwestern Atlantic Ocean and Carrib 


ean Sea 






Cape George, Nova Scotia 


1 


erosional 


Stea et al. (1998, 2001) 


Mark Clark, South Carolina 


1 


facies 


Cronin et al. (1981) 


Grape Bay, Bermuda 


2 


facies 


Muhs et al. (2002); Hearty et al. (2007) 


San Salvador Island, Bahamas 


3 


reef 


Chen et al. (1991) 


Great Inagua Island, Bahamas 


3 


reef; erosional 


Chen et al. (1991) 


Abaco Island, Bahamas 


3 


reef; erosional 


Hearty et al. (2007) 


Southern Barbados 


8 


reef 


Schellmann and Radtke (2004) 


Southwestern Atlantic Ocean 








Rio Grande do Sol coastal plain, Brazil 


1 


facies 


Tomazelli and Dillenburg (2007) 


Camarones, Patagonia, Argentina 


1 


erosional 


Rostami et al. (2000) 


Pacific Ocean 








Oahu, Hawaii 


3 


reef; corals; facies 


Hearty et al. (2007); Muhs et al. (2002) 


Mururoa Atoll 


1 


corals 


Camoin et al. (2001) 


Australia 








Eyre Peninsula 


1 


facies 


Murray- Wallace and Belperio (1991) 


Rottnest Island 


1 


reef 


Stirling et al. (1995); Hearty et al. (2007) 


Minim Cove 


1 


facies 


Hearty et al. (2007) 


Cape Range 


2 


reef 


Stirling et al. (1998) 


Houtman Abrohlos Islands 


8 


reef; facies; corals 


Zhu et al. (1993); Eisenhauer et al. (1996) 


Indian Ocean and Red Sea 








Red Sea 


30 


isotopic 


Rohling et al. (2008) 


KwaZulu-Natal, South Africa 


3 


erosional; facies 


Hobday (1975); Ramsay and Cooper (2002) 


Eastern Cape, South Africa 


1 


erosional 


Ramsay and Cooper (2002) 


Maldives Archipelago 


1 


facies 


Woodroffe (2005) 


La Digue Island, Seychelles 


2 


reef 


Israelson and Wohlfarth (1999) 


Aldabra Atoll, Seychelles 


3 


corals; facies 


Braithwaite et al. (1973) 


Polar regions 








Northern and Western Alaska 


3 


facies 


Brigham-Grette and Hopkins (1995) 


Wrangel Island, Siberia 


1 


facies 


Gualtieri et al. (2003) 


Western Spitsbergen 


3 


erosional 


Forman and Miller (1984); Andersson et al. (1999) 


Scoresby Sund, Greenland 


3 


facies 


Landvik et al. (1994); Vosgreau et al. (1994) 


Cape Ross, Antarctica 


1 


erosional 


Gardner et al. (2006) 
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by 350%, following the empirical observation of Scholz 
and Mangini (2007) of the overestimate of the preci- 
sion of ages from single sample measurements. When 
multiple measurements are reported, we employ their 
inverse- variance weighted mean. We expand the inverse- 
variance weighted standard deviation using a Student's 
t-distribution so that the 95% confidence interval spans 
±1.96(7, with (J the standard deviation, as in a Gaussian 
distribution. 

2.3. Tectonic uplift or thermal subsidence rate 

In order to remove the local tectonic contribution to 
paleo-sea level, we seek locally calibrated subsidence or 
uplift estimates for each locality. For most of the points 
in our database, no estimate of uplift or subsidence is 
available, but the value is expected to be near zero. For 
these location, we adopt an estimate for these locations 
of ± 1 cm/ky. In a few regions where estimates are 
available, including much of the Bahamas and Hawai'i, 
subsidence or uplift is on the order of 1-2 cm/ky. A 
few localities have exhibited uplift (Barbados, Patago- 
nia, southern England) or subsidence (the Netherlands, 
Pacific and Indian Ocean atolls) in excess of about 10 
cm/ky. The fastest uplifting locality in our database, 
Barbados, is rising at about 28 cm/ky. 




60 120 180 240 300 360 

Figure 1: Sites with at least one sea level indicator in our 
database. The symbol shapes reflect the nature of the indicators 
(upward triangles: isotopic; circles: reef terraces; downward tri- 
angles: coral biofacies; squares: sedimentary facies and non-coral 
biofacies; diamonds: erosional). The colors reflect the number of 
observations at a site (blue: 1; green: 2; magenta: 3; red: 4 or 
more). Locations marked by open symbols were excluded from 
our analysis because they have experienced strong near-fleld iso- 
static uplift for which our model cannot account. 



2.4- Coverage 

Our database attains fairly good geographic coverage, 
including the northwestern, northeastern, and south- 
western Atlantic coasts; the Carribean; Alaska, Green- 
land, Svalbard, and Siberia; Australia; the southwest- 
ern Indian coast; and Pacific and Indian Ocean islands 
(Figures 1 and 2; Table 1). (Because the physical model 
we employ for developing our covariance function does 
not account for the local effects of isostatic rebound, we 
were, however, unable to include near-field data from 
Greenland, Svalbard, and Antarctica in our analysis.) 
Where nearby localities subject to less uplift are avail- 
able, we have tried to limit the amount of data from 
rapidly uplifting sites, though we include Barbados be- 
cause of its prominence in the literature. However, given 
the long history of the geological study of Pleistocene 
sea level indicators , which began not long after the col- 
lapse of the Diluvian hypothesis in the early nineteenth 
century (e.g., Godwin- Austen, 1856), we do not claim 
that our database comprehensively represents the entire 
literature. 



Sea level database 

Position (r) 
Altitudes (z) 
Ages (t) 

Depositional range (D) 
Uplift/subsidence rate (u) 



Statistical 
Analysis 



Covariance Function from 
Allohistorical Physical Modeling 

Effects included: Eustasy, Gravitation, 
Isostatic mass compensation 

(Effects not yet included: elastic flexure, 
isostatic land motion, thermal expansion...) 



Initialization 

"True" sea level (f) 
"True" ages (g) 



I) Sea Level Measurements 

Calculation of measured sea level (s) 
* from z, D, u, g; imputation of censored 
data based on f and measurement 



3) Markov Chain Monte Carlo 
Sampling of Ages 

Stepwise sampling of new candidate 
"true" ages g 



2) Gaussian Process Regression 
for True Sea Level 

Regression of s for f using covariance 
function K(r,g) 



3. Statistical Model 

3.1. Preliminaries and Notation 

The ultimate goal of our statistical analysis is to 
determine the posterior probability distribution of sea 
level through time, conditioned upon the measurements 
in our database. Expressed symbolically, our aim is to 



Figure 3: Schematic illustration of the process used in our analy- 
sis. 
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error bars: — ±20 ky | ±10m 



Figure 2: Localities at which local sea level data exist in our database for time shces through the Last Interglacial. Points scale 
proportionately to the probability that they occur in the indicated interval. The horizontal lines are proportional to the standard 
deviation of the age measurement, and the vertical lines are proportional to the standard deviation of the sea level measurement . 
Censored data are indicated by upward (marine limiting) and downward (freshwater limiting) triangles. Color indicates the mean sea 
level estimate in meters. Open diamonds show near-field data points excluded from the analysis. 



evaluate the probability P(/(x, ^) |r, z, t, D, u) for loca- 
tions X on Earth's surface and times g, where / rep- 
resents the true value of sea level at x and g. In our 
database, each sea level indicator is assigned an index 
i = 1, . . . , and is characterized by 

r^, its exact geographic position, 

a noisy measurement of its altitude, 

ti^ a noisy measurement of its age, 

D^, a closed or open interval reflecting its deposi- 
tional range, and 

a noisy estimate of the long-term average uplift 
or subsidence rate. 

When is a closed interval, we replace it with d^, a 
Gaussian estimate of depositional depth characterized 
by the same mean and variance as the uniform distri- 
bution on D^, as discussed before. 



We collect these parameters into vectors r, z, t, D, 
u, and d. Similarly, we collect what will be the true sea 
levels in a vector f evaluated at the times g and locations 
X, whose elements /j, gj and Xj for j = 1, . . . , M are 
the desired sea levels and evaluation points. Only when 
geographical positions and depositional ranges are con- 
cerned does the bold vector notation serve double-duty: 
X and r are either coordinates or vectors of coordinates, 
and x^, Ti and Xj , Vj are individual sets of coordinates. 
Likewise, D is either a depositional range or an array 
of depositional ranges, and is an individual deposi- 
tional range. This dual purpose is not, however, likely 
to lead to confusion. 

3.2. Gaussian process regression 

We proceed from this point using a Gaussian pro- 
cess approach (Rasmussen and Williams, 2006). We 
must select some covariance function for true sea level, 
k{ri^ gi'^Vj^ gj), as we will address in section 3.4. Let 
(f , g) refer to the vectors of true sea levels and ages that 
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correspond to the vectors of measurements (z,t,D,u); 
i.e., with every entry we associate an entry 

{zi,ti^Y>i^Ui) for ah indices i = 1,...,A^. With the 
covariance function k given, we can then readily recover 
an estimate of true sea level at any arbitrary location 
and time g' through straight-forward kriging inter- 
polation (Press et al., 2007). We denote the mean and 
variance of this estimate by /(x',^') and V(/(x',^')), 
respectively. 

As before, the vectors f , x' and will collect the 
mean estimates of the sea levels at the desired points 
x' and in space and time. The sets of desired eval- 
uation points (x^-,^^), j = 1,...,M, and the measure- 
ments (vi^gi), ^ = 1, need not necessarily overlap. 
The matrix V collects the kriging (co) variance of f at 
and between (x',g'). Let K, K', and be the co- 
variances of (f,g) and/or (f^gO ^it the observed and 
desired points, i.e., let the symmetric square matrices 
K and K^' and the rectangular matrix be defined by 
their elements: 

Kij = k(vi,gi]Vj,gj) where i,j = l,...,iV, (1) 

Klj = fc(x-,^-;x^-,^^) where ij = 1,...,M, (2) 

Klj = k{Yi,gi]yij,g'j) where i = 1,...,N (3) 

and j 1,...,M. 

From this, the kriging step consists of calculating f , the 
M X 1 vector of mean sea level estimates at (x',g'), as 

f = K'"^K-^f , (4) 

which has 

V = K'' - K'^K-^K' (5) 

as its M X M covariance matrix. It is clear from the 
above that, when x' = r and g' = g, K = = K'^, 
and therefore f = f and = 0. In other words, when 
the queried points are identical to the measurement lo- 
cations, the interpolated values of true sea level remain 
unchanged and receive no kriging variance. 

We can therefore replace the problem of find- 
ing the posterior probability of sea level anywhere, 
P(/(x, ^)|r, z, t, D, u), with the more tractable problem 
of finding P(f , g|z, t, D, u), which is the posterior prob- 
ability of sea level at the smaller set of points defined by 
the measurement locations. After adjusting altitude Zi 
for uplift or subsidence rate Ui over a time gi, we define 
the corrected altitude z[ as 

z-=Zi-giUi, (6) 

with variance 
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and we define the sea level measurement Si and its vari- 
ance <j^^ as 

Si = zl-di, (8) 
^si=^l'i^(^di^ (9) 

where a^^, a^^, and a^- are the variances respectively of 
altitude Zi, uplift rate Ui, and depositional depth di. By 
B ayes' theorem, 

P(f,g|s,t)aP(s,t|f,g).P(f,g). (10) 

We drop the position variable r from the notation, since 
its values are fixed in the data set and implicit in the 
indexing of the other variables. For uncensored sea level 
measurements, we have the likelihood 

Pis,\fi,gi)^mi,'^i)- (11) 

In other words, the probability of observing sea level Si 
at a point in the data set that has a true sea level of fi is 
given by a Gaussian centered on the truth with variance 
cr^^. For censored data, 

Pisi\f.,9^)^mr,cri)-S{iz',-s,) e Di) (12) 

where S is an indicator function that is 1 when — Si 
is in the depositional range Di and otherwise. For 
instance, if is (— oo, —2], refiecting deposition at least 
two meters below mean tide level, then S would be 1 for 
Si ^ z • + 2 and otherwise. For age measurements, we 
have the likelihood 

P{U\g^)r.X{g,^al), (13) 

where cr|^ is the variance of age measurement ti . For the 
sea level vector f , we assume the prior 

P(f|g)~AA(0,K(g)), (14) 

where we use the notation K(g) for the covariance to 
emphasize its dependence on ages g. For the age vector 
g itself, we assume a uniform prior. 

3.3. Algorithm for sampling the sea level distribution 

To explore the distribution in equation 10, we use a 
recursive three-step algorithm (schematically illustrated 
in Figure 3) to generate updates of f , g, and s. We start 
by initializing g = t for all data points and = Zi—giUi 
and fi = Si = z[ — di for the uncensored ones. By simple 
kriging interpolation (equations 4 and 5), we estimate 
fi at the remaining data points. 

1. In step one of our algorithm, we calculate values 
of sea level measurements s from z, D, g and u. For 
uncensored data, Si is as defined in equation 8. For 
censored data, we sample Si from the distribution in 
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equation 12, with an additional variance term a'j-, the 
kriging variance of fi. 

2. In step two, we update our estimate of true sea level 
f based upon the new s as follows. We define the matrix 
of the sea level measurement noise N, with elements 
cr^^ along the diagonal and zero elsewhere. Then, by 
Gaussian process regression, paralleling equation 4, we 
calculate 



f = K(g)^(K(g) + N)-is, 



(15) 



the vector of sea level predictions and the vector of their 
variances 

S = diag{K(g)^(I - (K(g) + N)-iK(g))}, (16) 

where diag denotes the diagonal elements. 

3. In step three, we update our estimate of the true 
ages g. To do this, we follow a Markov Chain Monte 
Carlo approach applying the Metropolis-Hastings algo- 
rithm sequentially to each gi. Let represent g with 
element i removed. For each i, we sample from the 
distribution P(^i|t, g_i, f), which, by multiple applica- 
tions of Bayes' theorem and the facts that P(t|g) = 
UiPiUh) and that P(t|f) = P(t), reduces as 



P(^,|t,g_„f)(xP(t,|^0-^(f|g)-^(g). 



(17) 



The first term is given by equation 13, and the second 
term by equation 14. We can drop the third term be- 
cause of our assumption of a uniform prior for g. 

We generate test values using a Gaussian function 
q{g[]gi) centered at gi and bounded such that, when 
stratigraphic ordering is known, a point j that follows 
a point i always has gj ^ gi. (Where no bounds apply, 
q{a] h) = q{b; a).) For the sequences where relative ages 
are known more precisely than absolute ones, these are 
calculated in terms of time after the preceding point. 
Following the Metropolis-Hastings algorithm (Hastings, 
1970), we accept a candidate g[ with probability 



mm 



P(ff^|t,g_,,f)-g(g,;gO 
' P(5i|t,g_i,f) ■q{gi;gi) 

P{ti\gi)-P{i\g.u9i)-q{9'i;9i) 



min 1 



(18) 



So that we can assess results within a common temporal 
reference frame, we arbitrarily set the temporal variance 
cr|- for the first step of our longest quasi-continuous se- 
quence of data points (the sea level curve derived from 
the global oxygen isotope stack, for most runs) to zero. 

This algorithm, repeated a large number of times, 
samples the probability distribution described by equa- 
tion 10. We thin the results by storing every 20th sam- 
ple and account for burn-in by discarding the first 50 



stored samples. After several parallel executions of the 
algorithm, each of which store at least about 200 sam- 
ples, we check for convergence by inspecting the auto- 
correlation of stored values of g and discard executions 
that appear not to converge. To generate our target 
distribution P(/(x, ^)|s, r, t), we use kriging interpola- 
tion (equations 1-5) to estimate the sea level field at all 
spatial and temporal points of interest for each stored 
sample. 

We note that this algorithm, while satisfying from a 
theoretical perspective, could benefit from greater com- 
putational efficiency. The most time-consuming steps 
in its execution are the inversions of the covariance ma- 
trices, which for a database of n samples require 0{n^) 
operations. This inversion occurs once in step 2 and 
n + 1 times in step 3. Thus, each iteration of the algo- 
rithm is 0{n^). Repeating the algorithm a few thousand 
times in the courses of a Monte Carlo simulation with a 
database of about 150 points can therefore take several 
days; without increased efficiency, larger data sets will 
become unmanageable. 

3.4' The Covariance Function 

We use a covariance function that takes the form 

k{Yi,gi]Yj,gj) = ks{Yi]Yj) ■ h{gi]gj) (19) 

where kg and kt are respectively the spatial and tempo- 
ral components of the covariance function. To find suit- 
able ks and kt^ we employ a simple physical model of 
the gravitational effects of ice sheet melting and isostatic 
mass compensation. To generate /c^, we run the model 
300 times. During each run, we produce complete world 
maps of sea level at 20 random time slices, for a total of 
six thousand sea level maps. From these six thousand 
sea level maps, we compute the covariance among local 
sea levels at evenly spaced points on a Cartesian grid 
(spaced at 5° latitude and 10° longitude) and between 
these points, mean global sea level, and the volumes of 
the different ice sheets. We store the results as a lookup 
table, which is effectively over-sampled near the poles 
because of the use of the Cartesian grid. For the co- 
variance between two arbitrary points, we reference the 
closest grid points. 

Our physical model is based upon Woodward (1888) 
and the discussion thereof in Farrell and Clark (1976). 
Given a change in ice volume at a point p corresponding 
to global sea level rise of m^, the change in sea level Ri 
at an angular distance of is given by 



m 



\} 2sin(l9/2)} 



m, •(l+^(^)), 

1 I / PE 



(20) 
(21) 



where p^; is the mean density of the Earth (5.5 g/cm^) 
and is that of seawater (1.0 g/cm^). To prevent sin- 
gularities, we do not let R take values smaller than -10 
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m, a constraint equivalent to assuming that all points 
are at least 3° from a point mass. To approximate the 
gravitational effects of isostatic compensation, we place 
a slowly adjusting compensatory mass at p. The mass 
equivalent rric, in units of length, is given by the differ- 
ential equation 



d{mc — rrii) 



rur. 



dt 



(22) 



where Tc is the timescale of isostatic adjustment. This 
approach does not take into account the effects of flexure 
or isostasy on the solid Earth and so captures near-field 
behavior only crudely. In each model run, Tc is drawn 
from a uniform distribution between 0.1 and 14.1 ky. We 
include the low end of the distribution not because they 
are physically plausible isostatic timescales but because 
they allow us to include in the prior probability distri- 
bution scenarios in which eustasy completely dominates 
other processes. 

Including the gravitational effects of this compen- 
satory mass, the change in sea level at 6> is given 

by 



R{e) = mi + {nii - nic) ■ £(0). 



(23) 



We treat changes in sea level as resulting from changes 
in the mass of four ice sheets, representing the Lau- 
rentide (represented as two point ice sheets at 57.2° N, 
102.2°W and 56.5°N, 78.5°W), Greenland (at 65.5°N, 
49.5° W and 76.2°N, 22. 7° W), Scandinavian (at 64.2°N, 
14.5°E), and Antarctic (at 81.5°S, 176.5°W and 77.5°S, 
52. 5° W) ice sheets. These masses give rise to the grav- 
itational fingerprints shown in Figure 4. With long iso- 
static compensation timescales, these fingerprints are 
strongly expressed; with short isostatic compensation 
time scales, local sea level is nearly equal to global sea 
level. 

To generate reasonable ice sheet histories, we employ 
modifications of a modified form of the ICE-5G ice sheet 
history (Peltier, 2004) (Figure 5). In the base history, 
the ICE-5G ice sheet history from 21 ka to present is 
transposed to 146 to 125 ka, and occurs in reverse from 
125 to 104 ka. To generate a new scenario, the am- 
plitudes of melting associated with the Greenland, Eu- 
ropean, North American, and Antarctic ice sheets are 
each multiplied by uniformly distributed random num- 
bers between 0.6 and 1.4. For each ice sheet, the 21 
ky timescale of melting is stretched by a uniformly dis- 
tributed random factor drawn between 0.5 and 1.5 dur- 
ing the approach to the peak interglacial, and stretched 
by another random factor from the same distribution 
after the peak interglacial. 

In addition, we allow for up to 5 m of additional melt- 
ing from Antarctica between 127 and 123 ka, and up to 
7 m of melting from Greenland between 127 and 123 ka. 
The magnitudes are drawn from uniform random dis- 
tributions. As with the primary ice sheet melting, the 



time scales for the additional melting are stretched by 
random factors between 0.5 and 1.5. We also allow a 
brief lowstand at 125 ka caused by the loss of between 
0% and 100% of the additional melting from each ice 
sheet. 

While all sea levels are measured with respect to the 
modern sea level datum, changes in ice mass since the 
Last Glacial Maximum have not yet been fully compen- 
sated. In order for our model to report sea levels in 
this datum, it is therefore necessary to calculate the dif- 
ference between modern sea level and equilibrium sea 
level. For a given isostatic timescale Tc, we assume that 
ice masses were fully compensated (i.e., the mantle and 
ice sheets were in equilibrium) at 21 ka and follow the 
ICE-5G ice sheet histories to calculate modern sea level. 
In our random model runs, we assume that ice sheets 
are fully compensated at peak glaciation and subtract 
modern sea levels from all values. One portion of the 
resulting spatial covariance function, the covariance of 
local sea level with global sea level, is shown in Figure 
6. 

To find a suitable temporal covariance function kt, we 
explicitly calculate the temporal covariance function for 
global sea level and fit it using the sum of two Gaussian 
curves: 



9i - 9j 



(24) 



The best fitting values are r^^i = 4.31 ky, rk^2 = 0.41 
ky, /i = 0.97, and /s = 0.03 (Figure 11). 

Although our physical model is extremely simplified, 
neglecting all changes in the shape of the solid Earth 
as well as dynamic and rotational effects, and although 
coupling our statistical approach to a more sophisticated 
model may be useful, we do not think these simplifica- 
tions will significantly alter our global sea level projec- 
tions. We use our physical model to construct a rea- 
sonable prior probability distribution that discriminates 
regions where local sea level is strongly correlated with 
global sea level from regions where it is less strongly cor- 
related with global sea level. The simple model serves 
this purpose. Its failings may be more acute in evaluat- 
ing local sea level in the near-field of ice sheets, where 
the neglected changes to the shape of the solid Earth 
can be the dominant factor in sea level change. 

3.5. Validation of method using synthetic data 

To test our statistical model, we generated a synthetic 
sea level history in the same fashion as during the co- 
variance calibration process, using our base model with 
the arbitrary addition of 6.6 additional meters of melt 
from Greenland and 1.1 meters additional melt from 
Antarctica. We sampled the history at the same points 
in space and time (case A) as sampled by our database, 
under two conditions: (1) with no errors in ages and 10 
cm errors in sea level; (2) with the same chronological 
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and sea level errors as in the data set. We also ran con- 
ditions (1) and (2) excluding the oxygen-isotope derived 
global sea level curve (case B) and excluding both the 
oxygen-isotope curve and the Red Sea sea level curve of 
Rohling et al. (2008) (case C). We discuss the results of 
the validation analysis in section 4.1. 

3.6. Summary Statistics 

We report several summary statistics for each four- 
dimensional sea level distribution we discuss. We assess 
the median and quantiles of the data points by aggre- 
gating the sub-distributions determined by the stored 
mean and standard deviations of sea level data points 
and ages for each stored iteration of the MCMC model. 

We similarly assess the median and quantiles of global 
sea level by aggregating sub-distributions for global sea 
level over time from each stored iteration of the MCMC 
model. Each sub-distribution is determined from the 
stored mean and standard deviations of sea level data 
points and ages associated with the iteration by Gaus- 
sian process approximation of GSL at 500-year intervals 
from 115 to 140 ka. 

We calculate the 1000-year and 500-year average 
rates of global sea level change by taking the average 



slope of the global sea level curve from each iteration 
over 1000-year intervals or 500-year intervals and aggre- 
gating these curves to produce a distribution of global 
sea level rate. Note that these rates are average rates 
over several centuries; they place a lower bound on 
century-level rates of sea level rise. 

Of particular interest are the highest global sea level 
reached during the Last Interglacial and the fastest rate 
of sea level rise experienced. We report two sets of 
statistics relevant to these questions. First, we report 
the maximum of the median global sea level curve and 
sea level rate curve and its confidence interval. We 
also compute global sea level and global sea level rate 
exceedance probabilities. To do this, we sample each 
sub-distribution of global sea level one hundred times 
and aggregate all of these samples. In order to discount 
time points at which we have limited data, we incorpo- 
rate only the time points within each sample at which 
the posterior standard deviation is less than 15% of the 
prior standard deviation. We then identify the fraction 
of sample histories in which a 1000-year running average 
of GSL or the 1000-year or 500-year average rate exceeds 
a given value. The 95% probability exceedance levels. 
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ice volume scenarios used to build the covariance function k 




110 115 



Figure 5: Synthetic ice volume histories, expressed in terms of 
global sea level equivalent, used to build the covariance function. 
The histories are based on random distortions of the ICE-5G re- 
construction of LGM-to-modern ice volume. The pale curves il- 
lustrate the distribution of alternative histories used to generate 
the covariance function. The bold curves are used to generate the 
synthetic data for the validation analysis. 
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Figure 6: The spatial covariance of local sea level with global sea 
level ks(GSL;ri), normahzed to the variance of global sea level 
(c^GSL = 33.4 m). 



for instance, are the values that with 95% confidence 
we can say the sea level or sea level rate exceeded. For 
the rate maxima and rate exceedance probabilities, we 
focus on intervals beginning when global sea level was 
-10 m or higher, as we expect that ice sheet dynamics 
during these intervals will more closely resemble future 
ice sheet dynamics than will the behavior of ice sheets 
during intervals of lower GSL. 

We compute parallel summary statistics for Northern 
Hemisphere and Southern Hemisphere ice volume, ar- 
rived at by Gaussian process regression for these values 
instead of for GSL. 

To identify outliers among the data points, we com- 
pute the probability of a measurement given the assessed 
sea level distribution. To do this, we take the average 
over all N stored MCMC iterations of the probability 



that the parameter / (local sea level, global sea level, or 
age) with measured value fm ± cr^ was drawn from the 
distribution indicated by iteration i, with mean fi and 
standard deviation cji. For indicative points, the prob- 
ability for each iteration is given by a distribution 
with one degree of freedom on the parameter ^^^2^^2 • 
For limiting points, the probability is given by a cumula- 
tive normal distribution with mean fi — fm and variance 

To compare different probability distributions fi and 
f2 for a parameter / computed using different sub- 
sets of the data, we calculate the expected Mahalanobis 
distance. We sample with replacement 1000 pairs of 
MCMC iterations (fi,i,f2,i) from fi and f2. We then 
take the mean Mahalanobis distance between each pair. 
The Mahalanobis distance of each pair is given by 
V(fi,^ - f2,i)^(5:i,, + 5]2,^)(fi,^ - f2,i) where E,- , is the 
variance of fj^^. 



4. Results 

4.1. Validation analysis 

Comparisons of true and projected global sea level, 
rate of global sea level change, and Northern Hemi- 
sphere ice volume for the validation analysis are shown 
in Figure 7. (The cases with no sampling errors are 
shown in figure 12). Summary statistics are presented 
in Tables 4.1. We show both the 95% probability ex- 
ceedance levels and the exceedance probabilities corre- 
sponding to the true maximum values. The algorithm 
performs a good job of reconstructing global sea level, 
with the median projection often quite close to the true 
values, generally within the 67% confidence intervals, 
and always within the 95% confidence intervals. The 
same is true for rates and ice volumes. 

Cases Al, Bl, and CI highlight the interpretive limi- 
tations imposed by sampling (Tables 4.1 and Figure 12). 
All of these "perfect knowledge" cases do a good job of 
reconstructing global sea level and rates for the period 
of highest sea level. The maximum of the median pro- 
jections for GSL and rate are quite close to the true 
maxima, which fall between the 58% and 72% proba- 
bility exceedance levels. The true rates of change fall 
between the 23% and 64% probability exceedance val- 
ues. Without the d^^O curve included (as in cases Bl 
and CI), the resolution of the curves becomes poor be- 
fore about 130 ka, but in all cases the reconstructions 
do a good job of resolving details, including the 125 ka 
drawdown, after 130 ka. However, even under the best 
of circumstances, we do a mediocre job of reconstructing 
ice volumes; the 67% confidence intervals for our maxi- 
mum ice volume projections span about ±6 m. The true 
NH ice volumes fall between the 69% ans 72% probabil- 
ity exceedance levels, while the true SH ice volume falls 
between the 92% and 94% probability levels. 
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In cases A2, B2, and C2, where the measurement er- 
rors are included, they prevent resolution of some of the 
details of the curves, including the brief 125 ka draw- 
down (Tables 4.1 and Figure 7). Case B2 preserves a 
suggestion of the drawdown, while case A2 resolves re- 
solves multiple wobbles of sea level instead of a single 
drawdown. Case C2 finds a smoother sea level peak. 
Compared to cases Al-Cl, these cases exhibit increased 
uncertainty in estimates of the maximum of the median 
GSL and GSL rate curve, but preserve accuracy, with 
the true values remaining well within the 67% confi- 
dence interval. Compared to the more accurate cases, 
the exceedance levels are biased toward higher values. 
The true value of the GSL maximum falls between the 
79% and 91% probability exceedance levels, while the 
true value of the rate of change falls between the 83% 
and 91% probability exceedance values. The precision 
and accuracy of the ice volume projections are compara- 
ble to those of cases Al-Cl, indicating that insufficient 
sampling rather than measurement error is the major 
source of uncertainty in the ice projections. Based on 
the exceedance levels for these cases, we employ the 80% 
and 95% probability levels to bracket our estimates of 
maximum values in our analysis of the real data. 

4.2. Global analysis 

Applying our algorithm to the full data set of LIG sea 
level indicators (Tables 4.2 and Figure 8) reveals a peak 
median GSL of 4.8 =b 2.7 m (67% confidence interval) 
centered at 124 ka. The 95% and 80% probability ex- 
ceedance value for GSL are 5.8 and 6.8 m. This result 
is sensitive to the subset of the data examined (Figure 
13). Excluding the global S^^O curve but retaining the 
Red Sea curve yields a peak median GSL of 9.6 ±3.2 m, 
consistent with some of the high values that character- 
ize the Red Sea curve. (The Red Sea curve itself has a 
peak value of 12.4 ± 3.0 m (Icr).) The associated 95% 
and 80% probability exceedance values are 7.9 and 9.3 
m. Excluding both the oxygen isotope curve and the 
Red Sea curve yields a peak median GSL of 5.7 ±4.3 m 
and 95% and 80% probability exceedance values of 7.6 
and 9.1 m. We therefore conclude that global sea level 
during the Last Interglacial indeed reached significantly 
higher levels than present, probably in the range of 6-9 
m higher. 

The 95% and 80% probability exceedance values for 
1000-year average GSL rise rate during the interval 
when GSL was > —10 m range from 8.2 m/ky to 10.7 
m/ky. We emphasize that these values by no means ex- 
clude faster intervals of sea level rise lasting for a few 
centuries. 

As expected from our validation analysis, the data is 
insufficient to make strong statements about the source 
of the melt water that fed higher sea levels. They do, 
however, indicate that both the Northern Hemisphere 
ice sheets and the Antarctic ice sheets contributed. The 
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analysis of the full data set indicates that the minimum 
NH ice volume was 3.4±6.4 m GSL equivalent smaller 
than today and that the minimum SH ice volume was 
2.1 ± 5.7 m smaller than today. The ranges defined by 
the 95% and 80% probability exceedance values for each 
hemisphere are respectively 0.2-4.0 m and 0.3-3.5 m. 

The two subset analyses indicate smaller ice volumes, 
consistent with the higher global sea levels they also 
indicate. The case without the global oxygen isotope 
curve indicates a minimum NH ice volume of 7.5 ± 6.5 
m GSL equivalent and a minimum SH ice volume of 
2.0 ± 5.7 m GSL equivalent, with ranges defined by the 
95% and 80% probability exceedance values of 5.4-8.9 
m and 2.3-5.3 m. The case without both the global 
oxygen isotope curve and the Red Sea curve indicates 
a minimum NH ice volume of 4.5 ± 7.3 m but is oth- 
erwise nearly identical to the case with the Red Sea: a 
minimum SH ice volume of 1.8 ± 5.8 m, and ranges of 
of 5.5-9.0 m and 2.5-5.5 m. 

4.3. Cross-validation analysis of researcher bias 

One challenge in integrating local sea level records 
produced by many researchers is researcher bias. We 
therefore conducted a cross-validation study by classify- 
ing the data points by research group and repeating the 
global analysis on subsets of data in which we exclude 
all results from one research group. To compare the re- 
sults of analyses on different subsets, we calculated the 
Mahalanobis distances between them in terms of global 
sea level, as described in section 3.6. As controls, we cal- 
culated the distance between our full analysis and nine 
individual parallel runs of the Monte Carlo simulation 
on the full data set. 

When measured in terms of global sea level, only one 
subset was more different from the full analysis than all 
nine control runs: the subset excluding the global S^^O 
curve. Given the number of points in this curve, this 
result should be expected a priori. 

4.4' Outlier Analysis 

To search for outliers, we estimated the posterior 
probabilities for each of our sea level measurements and 
age measurements given the distribution at each point 
for sea level and age projected by our statistical model. 
We performed this analysis both upon the results pro- 
duced from the full data sets and the results produced 
from the subset excluding the 6^^0 curve and the Red 
Sea curve. 

No data point was a strong outlier, but three sites 
generated sea level measurement probabilities between 
0.15 and 0.33. Global sea level at 132, 131, and 129 ka 
was moderately higher than the values inferred from 
the global S^^O curve; inferred sea levels of — 58± 15 m, 
—35 ± 10 m and —40 ± 11 m (la) gave rise to projected 
GSL of -28 ± 15 m, -13±}i m and -27 ± 7 m. These 
estimates suggest a relationship between S^^O and sea 
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Table 2: Statistics from Validation Analysis 





Max. Global Sea Level 




GSL Rate at > -10 


m 






Median Projection 


95% Exceed. 


Prob. of 


Median Projection 


95% Exceed. 


Prob. of 


Case 


Time (ka) 


Max. (m) 


Level (m) 


True Value 


Time (ka) 


Max. (m/ky) 


Level (m/ky) 


True Value 


True 


127 or 123 


7.7 






124.5 


7.5 






Al 


123.5 


7.6 ± 0.6 


7.3 


58% 


124.5 






5.6 


64% 


Bl 


127 




6.4 


67% 


124.5 


7.3 ± 1.0 




4.6 


23% 


CI 


123 


8.2 ± 1.7 


6.3 


72% 


124.5 


6.6 ± 2.3 




5.0 


57% 


A2 


125.5 


7 7+2.4 


7.3 


91% 


136 


7.0±^;3 




7.0 


91% 


B2 


125 


8.6±i? 


7.0 


87% 


123.5 


5.8±4;8 




6.4 


83% 


C2 


125.5 




6.1 


79% 


127 


4.4±i;l 




6.9 


91% 




Max. Northern! Hemisphere Ice 




Max. Southern Hemisph 


ere Ice 




True 


6.6 






1.1 






Al 




6.0 ±5.9 


1.3 


69% 




2.3 ±5.6 




0.2 


92% 


Bl 




6.3 ±6.1 


1.6 


70% 




2.4 ±5.7 




0.6 


93% 


CI 




5.6 ± 6.1 


1.7 


72% 




2.7 ±5.7 




0.7 


94% 


A2 




5.6 ±6.6 


2.6 


78% 




1.7±5.8 




1.0 


95% 


B2 




6.7 ± 6.3 


2.2 


77% 




2.2 ± 5.8 




0.0 


92% 


C2 




4.0±^i 


2.5 


78% 




1.4 ± 6.0 




0.6 


93% 



67% confidence intervals shown. 



Table 3: Summary Statistics from Data Analysis 





Max. Global Sea Level 


GSL Rate at > -10 m 




Median Projection 


Exceed. Level (m) 


Median Projection 


Exceed. Level ( 


m/ky) 


Case 


Time (ka) Max. (m) 


95% 


80% 


Time (ka) Max. (m/ky) 


95% 




80% 


Full Data Set 


124.0 4.8 ±2.7 


5.8 


6.8 


126.0 6.2±^:^ 


9.0 




10.4 


No 6^^0 


123.5 9.6 ±3.2 


7.9 


9.3 


126.0 10. 3±^;^ 


9.2 




10.7 


No S^^O or RS 


129.5 5.7 ±4.3 


7.6 


9.1 


133.0 6.1±^;^ 


8.2 




10.1 




Max. Northern Hemisphere Ice 


Max. Southern Hemisphere Ice 


Full Data Set 


3.4 ±6.4 


0.2 


4.0 


2.1 ±5.7 


0.3 




3.5 


No 5^^0 


7.5 ±6.5 


5.4 


8.9 


2.0 ±5.7 


2.3 




5.3 


No 5^^0 or RS 


4.5 ± 7.3 


5.5 


9.0 


1.8 ± 5.8 


2.5 




5.5 



67% confidence intervals shown for the maximum of the median projections. 95% and 80% probability exceedance levels, which we employ as 
estimates of the maxima, are levels exceeded in 95% and 80% of all histories sampled from the estimated posterior distribution, respectively. 



13 



Sea Level During the Last Interglacial Kopp et al. 



case A2 case B2 case C2 




J i ' \ \ • h i ' \ • h r ' — \ V-- ' h 

110 120 130 140 110 120 130 140 110 120 130 140 

age (ka) age (ka) age (ka) 



Figure 7: Probability density plots (blue) of global sea level (top row), 1000-year average global sea level rates (middle row) and Northern 
Hemisphere ice volume (bottom row) as recovered by sampling a synthetic data set (black). The synthetic data was sampled with the 
same errors as in the real data set. Case A2 includes all data, case B2 excludes the global 6^^0 curve, and case C2 excludes both the 
global S^^O and the Red Sea data. The sohd hue marks the median projection, dashed lines mark the 16th and 84th percentiles, and 
dotted lines mark the 2.5th and 97.5th percentiles. Crosses at bottom mark the median posterior estimates of the sample ages. 



level of about -30 m/%o, a somewhat surprising result 
because such low slopes are more characteristic of glacia- 
tions than deglaciations (Waelbroeck et al., 2002). 

At Kahe Beach State Park, Oahu, Hawai'i, (Hearty 
et al., 2007) describe a marine conglomerate at 12 m 
above present sea level. Corrected for uplift of Oahu, 
this suggests a paleo-sea level of at least 9.6 ± 1.3(lcr) 
m. Our model instead assigns a sea level of 7.5 ± 1.7 m. 
When the S^^O and Red Sea curves are not included, it 
assigns a sea level of 8.0±2;4 These results suggest 
that the sea level indicator should be reexamined. 

Finally, our model identifies as an outlier early We- 
ichselian (post-Eemian) lacustrine sediment from a bor- 
ing in the North Sea (Zagwijn, 1983). The sediment 
indicates freshwater conditions at a relative sea level 
of about -40 m, which we adjust to — 23 ± 3 m based 
upon the subsidence estimates of Kooi et al. (1998). The 
model, however, places sea level at — 18.5±8;9 m. In the 
absence of the S^^O and Red Sea curves, it places sea 



level at —10.4 ± 4.4 m. In the absence of these curves, 
the model also identifies as an outlier early Eemian la- 
custrine sediment in the same core. These results sug- 
gest that the North Sea in the region of this boring is 
subsiding faster than the Kooi et al. (1998) estimates. 

5. Discussion 

5.1. Comparison to past estimates of Last Interglacial 
sea level 

The Fourth Assessment Report of the IPCC (Jansen 
et al., 2007) estimates that LIG sea level reached values 
of 4-6 m above present. Their estimate is qualitatively 
similar to the result of our analysis, but ascribes an 
excessive degree of precision and may be somewhat low. 
Our analysis suggest that sea level peaked during the 
Last Interglacial between about 5 and 9 meters above 
present. A more thorough analysis of sea level records. 
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Figure 8: Probability density plots of global sea level, 1000-year average global sea level rates, Northern Hemisphere ice volume and 
Southern Hemisphere ice volume projected from our fuh database. Dashed hues mark the 16th and 84th percentiles; dotted lines mark 
the 2.5th and 97.5th percentiles. Crosses at bottom mark the median posterior estimates of the sample ages. 



in other words, seems to permit even higher sea level 
than the IPCC's figure. 

The IPCC summary estimate, like other similar esti- 
mates (e.g., Overpeck et al., 2006), is based upon ex- 
amining a few key coral reef terrace localities. The 
IPCC highlights Hawaii and Bermuda (Muhs et al., 
2002); Overpeck et al. (2006) also highlight the Ba- 
hamas, Western Australia, and the Seychelles Islands. 
All these localities are basically tectonically stable and 
experience slow thermal subsidence. If one had to draw 
conclusions about global sea level from a small number 
of local sea level measurements, these would be reason- 
able sites at which to look. 

Other commonly considered localities, such as Barba- 
dos (e.g., Schellmann and Radtke, 2004) and the Huon 
Peninsula (Esat et al., 1999), are rapidly uplifting lo- 
calities. These sites have advantages as relative sea 
level recorders, most notably that terraces recording sea 
levels below present are readily accessible. Assuming 
these sites have experienced a steady rate of uplift, they 
can help uncover sea level variations over fairly short 
timescales. However, they are poor sites from which to 
draw conclusions about absolute sea levels, as recover- 



ing this information requires a precise estimate of uplift 
rate. 



To our knowledge, no previous author has accounted 
for the effects of glacial isostatic adjustment (GIA) in 
drawing conclusions about global sea level and ice vol- 
ume from Last Interglacial sea level records. As Lam- 
beck and Nakada (1992) demonstrated, understanding 
the infiuence of these effects is critical. Without this un- 
derstanding, local sea level highstands could be falsely 
interpreted as refiecting global highstands. In the mid- 
to-late Holocene, for instance, GIA has produced local 
highstands in far-field equatorial islands of about 1-3 
m above present levels (Mitrovica and Peltier, 1991); 
looking only at these sites in isolation, one might falsely 
infer that global sea level was higher and global ice vol- 
ume significantly smaller in the mid-Holocene than at 
present. Our statistical model uses the covariance be- 
tween local and global sea level to correct for these ef- 
fects; our results indicate that the apparent high Last 
Interglacial global sea levels are real. 
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5.2. ''Fingerprinting^^ analysis of meltwater sources 

Just as our model can predict global sea level from 
local sea level measurements, it can also predict changes 
in ice sheet volumes. However, as the validation analysis 
showed, the current database combined with the simple 
GIA model we employ cannot predict these changes with 
great precision (Figure 8c, d). 

The analysis does suggest that it is likely both that 
melting of northern hemisphere ice sheets and also melt- 
ing of the Antarctic ice sheet contributed to higher sea 
levels. More near-field sea level measurements would 
help increase the precision of these estimates; however, 
incorporating near-field measurements into the analy- 
sis requires a more sophisticated physical model. Our 
model correctly accounts for neither isostatic uplift nor 
fiexure, both of which are important in interpreting 
near-field sites. Indeed, even though we include some 
near-field data from Svalbard and Greenland in our 
database, we could not incorporate these points into 
our analysis. 

However, our approach could readily be adapted to 
incorporate a sophisticated glacial isostatic adjustment 
model like that of Mitrovica and Milne (2003). Such a 
model could be used to generate a spatial and temporal 
covariance function in place of our simple model. To 
do so would require a method for randomly generating 
hundreds of plausible ice sheet histories extending from 
before the Last Interglacial to today and modeling the 
full gravitational, elastic, and isostatic effects of these 
changing ice sheet volumes on local sea level around the 
globe. Such an effort would be computational much 
more intensive than employing our simple model, which 
takes less than fifteen minutes on a desktop computer 
to run 300 simulations. Nonetheless, it is quite feasible. 

5.3. The need for more data and opportunities for future 
research 

Mapping the uncertainty of local sea level projections 
highlights the regions whence more data is most strongly 
needed. The "Data Need Index" (Figure 9) is the mean 
of the ratio of the posterior standard deviation to the 
prior standard deviation over the time period between 
118 and 130 ka. The greatest need is in the near-field 
of the major ice sheets; as noted above, the need in 
these regions is to both expand the database and em- 
ploy a more sophisticated GIA model in our algorithm. 
The second greatest need for data is in the southern 
Indian Ocean, the regional antipodal to the Laurentide 
Ice Sheet, and thus more sensitive than the rest of the 
intermediate and far-field to GIA effects associated with 
Laurentide melting. Studies in areas such as Madagas- 
car, Mauritius, and the islands of the French Southern 
and Antarctic Lands would be helpful in this regard. 

In compiling the LIG sea level database, we also found 
a number of region where sea level indicators require 
further investigation. For instance, although Britain is 
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Figure 9: Map of the data need index. We calculate this index by 
averaging the ratio of the posterior standard deviation to the prior 
standard deviation over the time period between 118 and 130 ka. 
Units are percentages. Dots indicate the positions of observations. 
Open dots are observations not included in the analysis. 



on a tectonically stable passive margin, erosional ter- 
races appear to get progressively older with increasing 
elevation. Westaway et al. (2006) estimated Pleistocene 
uplift rates in the vicinity of the Solent river system 
range of ~ 10 m/ky. The causes of this uplift are un- 
certain, but might be linked to isostatic effects caused 
by erosional unroofing and the transport of sediment 
from continent to slope. A simple isostatic calculation 
indicates this method requires the removal of ~50 m of 
sediment per 100 ky. Clayton (1996) estimates that an 
average thickness of ~145 m of sediment was removed 
from the land of the British Isles to the continental shelf 
during the last glaciation; this removal could therefore 
be a potential cause. Because the British Isles are in a 
crucial region to look for the fingerprint of Greenland 
melting, a better understanding of regional uplift would 
be extremely helpful 

Braithwaite (1984) described numerous terraces in the 
coastal limestone of Kenya which range in elevation from 
-35 m to +20 m but lack good age constraints. These 
represent ready targets for modern dating techniques. 

5.4' Rates of sea level change 

we expect that ice sheet dynamics during these inter- 
vals will more closely resemble future ice sheet dynamics 
than will the behavior of ice sheets during intervals of 
lower GSL. 

Our results suggest that, during the interval of the 
Last Interglacial when sea level was — 10 m or higher, the 
rate of sea level rise, averaged over one thousand years, 
reached values of at least about 8-11 m/ky. Within the 
resolution of our data, we cannot confidently resolve 
rates of sea level change over shorter periods of time. 
These values are consistent with Carlson et al. (2008) 's 
estimates of the rate of the contribution of the Lauren- 
tide Ice Sheet meltwater to global sea level during the 
early Holocene; they estimate a Laurentide contribution 
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of about 7 m/ky during the period when global sea level 
climbed above —10 m. 

Ice volume during the late deglacial rise at the start 
of the Last Interglacial was only slightly larger than at 
present. The Laurentide Ice Sheet would have been 
a shrunken remnant of its once extensive mass - or, 
perhaps two small remanents, one over Quebec and 
Labrador and one over eastern Nunavut and Baffin Is- 
land, as in the late early Holocene (Carlson et al., 2008). 
It was within a factor of two in size of the present Green- 
land Ice Sheet, and its dynamics may therefore have 
been analogous to those of the Greenland Ice Sheet. 
Given a sufficient forcing, the results from the Last 
Interglacial suggest that the present ice sheets could 
sustain a rate of global sea level rise of about 80-110 
cm/century for several centuries, with these rates po- 
tentially spiking to higher values for shorter periods. 

6. Conclusion 

Contrary to an analogy commonly taught in introduc- 
tory classes, adding water from melting land ice to the 
ocean is not like pouring water into a bathtub. Many 
factors other than the changing volume of water in the 
ocean modulate the influence of melting ice sheets on 
local sea level. These factor include: the effects of the 
distribution of ice, water, and sediment on the geoid; 
lithospheric flexure; isostatic adjustment; and tectonic 
uplift and subsidence, as well as dynamic effects, which 
are of lesser concern on multi-century timescales. 

Consequently, global sea level and global ice volume 
cannot be accurately inferred simply by examining local 
sea level at one or two localities, yet this is the path most 
commonly taken when discussing the Last Interglacial. 
Our approach, which combines an extensive database 
with a new statistical algorithm for analyzing quanti- 
tative paleoenvironmental data with both interpretive 
and geochronological errors, offers better control. The 
results of our analysis support the common hypothesis 
that Last Interglacial global sea level was higher than 
present. We flnd that peak GSL was probably 6 to 9 m 
higher than present. 

The Last Interglacial was only slightly warmer than 
present, with polar temperatures similar to those ex- 
pected under a low-end, ~ 2°C warming scenario. 
Nonetheless, it appears to have been associated with 
substantially smaller ice sheets than exist at present. 
Our results indicate that both the Greenland Ice Sheet 
and the Antarctic Ice Sheets were at least smaller than 
today. Global sea levels flve meters higher than present 
could have been produced by a signiflcantly smaller 
Greenland Ice Sheet, a signiflcantly smaller Antarctic 
Ice Sheet, or both. Global sea levels nine meters higher 
than present would have required signiflcantly smaller 
ice sheets in both hemisphere, including nearly complete 
melting of either the Greenland and the West Antarc- 
tic Ice Sheets. This paleoclimatic constraint emphasizes 
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the vulnerability of ice sheets to even relatively low lev- 
els of sustained global warming. 
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A. Database 

Our database of Last Interglacial sea level indicators 
is based on a literature search for indicators with best 
estimates of ages between 140 and 110 ka. We charac- 
terized each candidate indicator by several quantitative 
parameters: 

• Location (latitude and longitude) 

• Age estimate 

• Altitude of indicator 

• Indicative meaning (i.e., the relationship between 
indicator and sea level) 

• Estimate of tectonic uplift or subsidence rate 

• Stratigraphic order, when more than one point 
comes from the same site; where available, we also 
include estimates of the relative ages of points at 
the same site 

The database is recorded in a spreadsheet that accom- 
panies this supplemental material. Three of the sites are 
re-analyses of data available elsewhere that require spe- 
cial explication: the sea level curve derived from the 
Lisiecki and Raymo (2005) global oxygen isotope curve, 
the re- aligned Red Sea sea level curve of Rohling et al. 
(2008), and a subsidence-corrected Dutch sea level curve 
based on Zagwijn (1983). 
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A.l. Global oxygen isotope stack 

By making assumptions about the isotopic composi- 
tion of ice sheets and a paleotemperature constraint, it 
is possible to estimate ice volume and thus global sea 
level from carbonate S^^O. These records most com- 
monly come from the shells of foraminifera, as is the 
case for the global oxygen isotope stack of Lisiecki 
and Raymo (2005) (LR04). Independent paleotemper- 
ature estimates can be derived from the Mg/Ca ratio 
of foraminifera (e.g., Lea et al., 2002), alkenone com- 
position (e.g., Haywood et al., 2005), or clumped iso- 
topes (Ghosh et al., 2006). Alternatively, a sea level 
record can be derived by modeling or making assump- 
tions about the relationship between ice sheet volume, 
changes in ice sheet isotopic composition, and temper- 
ature (e.g., Bintanja et al., 2005). 

The simplest way to convert these oxygen isotopic 
measurements into a measure of sea level is to calibrate 
the curve against modern and Last Glacial Maximum 
values. By definition, modern global sea level is zero me- 
ters; with reference to the VPDB oxygen isotope stan- 
dard, the modern benthic foraminifera oxygen isotopic 
composition from the LR04 stack is 3.23 ± 0.06%o. The 
Last Glacial Maximum value (at 18 ka) is 5.08±0.06%o, 
while sea level was ~ —125 m (Peltier, 2004, e.g.,). 
From this result, one would derive a conversion fac- 
tor of ~ —67.6 ± 7.2 m/%0. Based on comparison 
of coral records with oxygen isotopes, however, Wael- 
broeck et al. (2002) observe that this conversion is not 
always linear. During the onset of glaciations, oxygen 
isotopes tend to respond more quickly than sea level, 
while during terminal deglaciations, sea level responds 
more quickly. These changes refiect changes in both ice 
sheet composition and the relationship between temper- 
ature and ice sheet melt. Near modern sea levels. North 
Atlantic records indicate a relationship during glacia- 
tion of ~ 30 m/%0, while during glacial terminations, 
the relationship is ~ 90 m/%o. We therefore use a re- 
lationship of 60 ± 30 m/%0 in our calculations. When 
this uncertainty is combined with measurement uncer- 
tainty, it yield a peak Last Interglacial value of 7.8 ±8.0 
m (from an oxygen isotopic composition at 123 ka of 
3.10±0.10%o. 

The absolute ages of the LR04 stack between 22 and 
120 ka are derived by alignment against the oxygen iso- 
topic record of the MD95-2042 core, from the Iberian 
margin (Shackleton et al., 2000). This core is, in turn, 
aligned between and 67 ka against the oxygen iso- 
topic record of the GRIP ice core (Johnsen et al., 1992). 
Termination II is assigned to start after ~135 ka based 
upon U-Th dating of corals terraces from Papua New 
Guinea (Stein et al., 1993). While this age model is 
not necessarily superior to alternative age models (for 
instance, that employed by Rohling et al. (2008)), we 
have aligned the other quasi-continuous records against 
it so as to provide a common reference frame. 
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Figure 10: The LR04 global sea level curve and local sea level 
curves from the Red Sea and the Netherlands. Dashed lines show 
la confidence intervals in sea level. The initial best alignment of 
the three curves is shown. 

A. 2. Red Sea 

The Red Sea record is a planktonic foraminiferal oxy- 
gen isotope record that takes advantage of the particular 
hydrology of the Red Sea (Siddall et al., 2003) and is 
therefore essentially a local record of sea level at the 
Strait of Bab el Mandab. The oxygen isotopic composi- 
tion of Red Sea water is controlled primarily by evapora- 
tion. Water exchange occurs between the Red Sea and 
the Indian Ocean occurs through the strait; when sea 
level is lower, water exchange decreases, which increases 
the residence time of water in the Red Sea and thus 
yields heavier oxygen isotope values. This greatly mag- 
nifies the isotopic effects of sea level change. The dif- 
ference between the modern and the Last Glacial Max- 
imum in the Red Sea is nearly 6%o, whereas in the open 
ocean the difference is approximately 1.8%o. 

Using a hydrological model, Rohling et al. (2008) con- 
structed a sea level record with a raw la precision of 
6 m for the Last Interglacial from two Red Sea cores 
sampled for oxygen isotopes at 10 cm resolution. They 
aligned their record temporally with the record derived 
from U/Th-dated Barbados coral data (Thompson and 
Goldstein, 2005); in this age model, their record has a 
temporal resolution of 200-400 years. It indicates that 
local sea level rose to at least 6 ± 3.5 m, and perhaps 
as high as 11 m, during the peak interglacial. 

We have for consistency realigned the Rohling et al. 
(2008) against the age model for the global oxygen iso- 
tope stack of Lisiecki and Raymo (2005), which is based 
primarily on alignment against the GRIP ice core. This 
realignment required shifting the curve earlier by 2.4 ka 
and expanding the duration between measurements 1.2 
times. We include in our database the re-aligned sea 
level curve derived from the KLll core, while Rohling 
et al. (2008) provide a higher resolution record than the 
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Figure 11: The temporal covariance function of global sea level 
(black) and the Gaussian fit used to approximate it in our model 
(red). 

KL09 core. 

A. 3. Netherlands 

The Dutch Eemian sea level record of Zagwijn (1983) 
is based on sedimentological and micropaleontological 
data from numerous cores through the Amsterdam and 
Amersfoort basins, as well as cores along the Noord- 
Holland coast, in Friesland, and in the North Sea. Sea 
level indicators in these cores are provided by facies 
transitions representing, for example, the infiltration of 
marine water into a freshwater lake or the maximum el- 
evation of clays deposited in a salt-marsh environment. 
Relative age constraints are provided by characteristic 
Eemian pollen zones, many of which have durations es- 
tablished to fairly high precision based upon the count- 
ing of varves in an annually-layered lacustrine diatomite 
in northwestern Germany (Zagwijn, 1996). We place 
peak sea level in the middle third of zone E5 based 
upon the position of the maximum fiooding interval 
within the more recent Amsterdam- Terminal borehole 
(van Leeuwen et al., 2000). We estimate absolute ages 
from these relative ages by aligning the sea level curve 
against the Lisiecki and Raymo (2005) global oxygen 
isotope stack. 

Zagwijn reported sea level estimates without correc- 
tion for long-term isostasy, compaction, or tectonics. 
To correct for these factors, we use the backst ripping- 
derived Quaternary rate estimates of Kooi et al. (1998). 
These vary considerably across the Netherlands and the 
North Sea, ranging from about 12 cm/ky in Amersfoort 
to about 18 cm/ky in Petten. Thus adjusted, Zagwijn's 
data indicate that a maximum local sea level of about 
5 ± 2 m was attained in the Netherlands for much of the 
Last Interglacial. 
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Figure 12: Probability density plots (blue) of global sea level (top row), 1000-year average global sea level rates (middle row) and 
Northern Hemisphere ice volume (bottom row) as recovered by sampling a synthetic data set (black). The synthetic data was sampled 
with no temporal errors and 10 cm elevation errors. Case Al includes ah data, case Bl excludes the global S^^O curve, and case CI 
excludes both the global 6^^0 and the Red Sea data. Dashed lines mark the 16th and 84th percentiles; dotted lines mark the 2.5th and 
97.5th percentiles. Crosses at bottom mark the median posterior estimates of the sample ages. 
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Figure 13: Probability density plots (blue) of global sea level (top row), 1000-year average global sea level rates (middle row) and 
Northern Hemisphere ice volume (bottom row) projected from our full database (left column), our database excluding the global 6^^0, 
and our database excluding both the global S^^O and the Red Sea data. Dashed hues mark the 16th and 84th percentiles; dotted Hues 
mark the 2.5th and 97.5th percentiles. Crosses at bottom mark the median posterior estimates of the sample ages. 
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